function frame_data = h5load(h5file, vargin)
temp = hdf5info(h5file);
x = temp.GroupHierarchy;
frame_data = struct();

if nargin == 2
    Index = vargin;
    dataSize = min( length(x.Groups), length(Index) );
else
    dataSize = length(x.Groups);
    Index = 1:dataSize;
end

% all the simulation step
%for i =  1:length(x.Groups)     % loop over all the simulation steps
for iter=1:dataSize
    i=Index(iter);
    frame = x.Groups(i);
    frame_name = sprintf('frame%06d', i);
    %frame_name = x.Groups(i).Name(1:end);
    frame_data.(frame_name) = struct();
    
    % Dataset  stepsize on the same level of contacts & constraints
    for n = 1 : length(x.Groups(i).Datasets)
        % Hierarchy name with the absolute path
        DatasetHirName  = x.Groups(i).Datasets(n).Name;
        % pure name to save, here is specially for the step_size
        DatasetPureName = DatasetHirName(length(x.Groups(i).Name)+2 : end);
        if strcmp(DatasetPureName, 'timestep')
            DatasetPureName = 'step_size';
        end
        frame_data.(frame_name).(DatasetPureName) = h5read(x.Filename, DatasetHirName);
    end
    % each step with bodies, constraints and contacts
    for j =  1:length(frame.Groups)  %
        group = frame.Groups(j);
        g_name = group.Name(length(frame.Name)+2:end);
        frame_data.(frame_name).(g_name) = struct();
        % each field as datasest on the lowest level
        for k =  1:length(group.Datasets)
            data = group.Datasets(k);
            d_name = data.Name(length(group.Name)+2 : end);
            f = h5read(x.Filename, data.Name);
            f = f';
            [m, n]  = size(f);
            if m==1 && m < n
                f = f';
            end
            % contacts: all the PSI are saved as: POSITIVE for no
            % penetration NEGATIVE for penetration
         
            if strcmp(d_name, 'gaps')
                %the case when penetration depth is saved as positive
                d_name = 'PSI';
            elseif strcmp(d_name, 'world_inertia_tensor')
                d_name = 'inertia';
            elseif strcmp(d_name, 'h') || strcmp(d_name, 'timestep')
                d_name = 'step_size';
            end
            frame_data.(frame_name).(g_name).(d_name) = f;
        end % end for k
    end % end for
end  % end for iter
end

